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This  memorandum  describes  NUSL  Program  S1442  written  in  FORTRAN  V. 
This  program  uses  normal  mode  theory > to  predict  acoustic  propagation 
over  a flat  homogeneous  ocean  bottom. 


Program  SI 441  is  an  extension  of  a program  written  for  NUSL  by 
A.D.  Little  Inc.  (Reference  1 and  2).  In  the  A.D.  Little  program  the 
amplitude  distribution  of  an  acoustic  signal  as  a function  of  depth  is 
determined  for  a given  mode  by  means  of  the  numerical  solution  of  the 
acoustic  wave  equation  for  given  boundary  conditions.  Program  S1441 
calculates  and  produces  Calcomp  plots  of  the  following  for  any  mode, 
frequency  and  velocity  profile: 


a. 


any  mode. 


Amplitude  as  a function  of  depth  and  the  ray  equivalent  of 


b.  Group  velocity,  phase  velocity,  excitation  pressure,  and 
angle  of  incidence  of  sound  waves  striking  the  boundaries. 


c. 


Propagation  loss  as  a function  of  range. 
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Normal  mode  theory  la  based  on  the  assumption  that  at  large  distances 
from  a source  standing  waves  formed  by  energy  striking  the  boundaries  In 
certain  preferred  directions  form  the  main  part  of  the  sound  field. 

Solution  of  Wave  Equation 

There  are  three  methods  by  which  these  preferred  directions  and  the 
amplitude  distribution  of  the  standing  waves  may  be  determined.  First 
one  may  find  the  solution  in  closed  form  of  the  wave  equation. 

(1> 

Where  Is  a displacement  or  velocity  potential  and  C is  the 

velocity  of  sound  in  the  medium  considered.  This  is  done  in  Reference 
3 by  integrating  equation  (l)  subject  to  the  boundary  conditions  in  the 
complex  plane  and  approximating  the  normal  mode  solution  for  the  standing 
waves  by  the  evaluation  obtained  from  the  residues  in  the  integration. 

A second  method  is  the  solution  of  equation  (1)  by  direct  numerical 
integration.  This  was  done  by  A. D.  Little  Inc  by  means  of  a computer 
program  described  in  References  1 and  2. 

Third  it  has  been  shown  in  Reference  4 that  one  may  consider  the  ray 
paths,  which  undergo  total  reflection  at  the  boundaries  and  whose 
successive  upgoing  and  downing  rays  are  in  phase,  as  the  descriptions  of 
the  modes. 


A.  Numerical  Integration  of  the  Wave  Equation 

The  basis  of  this  program  is  the  direct  numerical  integration  of  the 
wave  equation.  This  method  which  is  extremely  well  suited  for  analysis 
by  means  of  a high  speed  computer  is  discussed  below. 

When  equation  (1)  is  solved  for  C^>  the  incremental  pressure  p is 
found  by  definition 

p ’ {2) 
is  the  displacement  potential 

is  the  material  density  at  the  point  considered 


Where  <|> 

f 
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It  is  found 
depth  l so  that 


that  equation  (1)  is  separable  in  terms  of  range  r and 
the  pressure  amplitude  can  be  written  as 


p.N  s FCr  V uCl'J 


(3) 


and  a differential 


4^ 


equation  in  terms  of 


where  vO  - the  radial  frequency 


may  be  obtained. 


(4) 


c = sound  velocity 

u = u(z)  is  the  pressure  amplitude  distribution  as  a function 
of  depth 

horizontal  wave  number 

Vsr  » ^ !»vv\  <S>  (5) 

where  © is  measured  relative  to  the  normal  to  the  ocean  bottom. 

The  physical  picture  presented  by  equation  (3)  is  that  of  a 
standing  wave  with  a particular  pressure  amplitude  distribution  as  a 
function  of  depth, mC*),  which  travels  unchanged  in  shape  as  it  progresses 
in  the  r direction. 


Equution  (4)  may  be 


written  as 

« 


o 


where 


(6) 

(7) 


In  solving  equation  (4)  there  are  three  constraints  which  must  be 
imposed  upon  any  solution.  These  constraints  are  due  to  boundary 
conditions  at  the  interface  between  the  water  and  bottom  material  and  at 
the  waters  surface. 


First  due  to  the  necessity  of  continuity  of  pressure  and  particle 
velocity  in  the  vertical  direction  at  the  bottom  interface,  it  is  necessary 
that  J 


\ du> 

'“'i  c\  1 


(Sa) 
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where  the  terms  having  subscript  b refer  to  these  quantities  in  the 
bottom  material  on  one  side  of  the  interface  while  those  with  subscript 
1 refer  to  the  water  side  of  the  interface. 


Second,  since  we  assume  a pressure  release  surface  at  the  air-water 
boundary  then 


U»  O 


(8b) 


at  this  surface. 

Third,  the  modes  by  definition  involve  propagation  with  energy  which 
strikes  the  bottom  at  angles  larger  than  the  critical  angle  so  that 
"total  reflection"  occurs.  Since  the  critical  angle  measured  relative 
to  the  normal  to  the  interface  is  given  by  %,'nw  0 » 

° Ov> 

(naturally  cto  must  be  creater  thane,) 
then 

to 

must  hold  for  all  energy  striking  the  bottom. 

Since 

Vr  5 ^ Si\v»  & 

using  euation  (7)  and  the  condition  of  continuity  of  pressure  we  obtain 
the  result  that 

\ ^ ^ (gc) 

C* 

at  the  bottom  interface. 

Upon  examining  equation  (6)  it  can  be  seen  that  a solution  of  the 
equation  is  of  the  form  

- oe 

for  given  values  of  Sa'-J  and  k,-  • B is  a constant.  If$(i)is 
positive  Uj)  is  in  the  form  of  a sinosoid.  This  kind  of  solution  is 
obtained  for  the  interference  pattern  between  upgoing  and  downgoing 
waves  in  the  water.  If  j(A)  is  negative  uii)  is  in  the  form  of  an 
exponential.  is  negative  for  the  following  two  cases.  It  is 

negative  everywhere  in  the  bottom  and  it  is  negative  in  the  water  at 
depths  which  correspond  to  shadow  zones  caused  by  vertexing  of  rays 
which  form  a mode.  Both  distributions  are  a result  of  the  condition  of 
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continuity  of  pressure  In  the  medium.  Thus  when  there  is  "total 
reflection"  at  a level  either  by  vertexing  or  reflection  from  a boundary 
and  the  pressure  is  finite  at  that  level  then  at  adjacent  levels  there 
may  be  a decay  (whose  rate  is  determined  by  the  boundary  conditions), 
but  not  a discontinuous  step  to  zero  pressure.  An  example  of  a distribu- 
tion Involving  both  exponentials  is  shown  in  Figure  1 which  shows  the 
amplitude  distribution  as  a function  of  depth  and  Figure  2 which  shows 
the  ray  equivalent,  of  the  mode.  In  Figure  2 it  is  seen  that  the  ray 
equivalent  of  the  normal  mode  vertexes  at  a depth  of  43  feet.  Therefore 
In  Figure  1 the  pressure  amplitude  distribution  is  in  the  form  of 
an  exponential  between  the  surface  of  the  ocean  and  a depth  of  43  feet. 

It  can  also  be  seen  in  Figure  1 that  the  pressure  amplitude  distribution 
in  the  bottom  is  in  the  form  of  an  exponential. 


In  order  to  Integrate  equation  (6)  numerically  formulas  relating 
'A*,,  and  its  derive tivevX»,x with  the  quantities  U*  and  must  be  obtained 
(The  subscripts  signify  the  depth  at  which  they  are  calculated).  This 


and 


Is  done  by  writing  a Taylor  series  for 
this  procedure  is  given  in  References  1 and  2.  The  results  are 
summarized  in  equations  (9)  and  (10) 

* - \ 


The  details  of 


(9) 


u 


A*  \ 


Vi 


u 


*{10) 


lv\.  y 


where 

V*  is  the  increment  between  level  v\  and  level  v\y \ . 
and  ivyv>  are  the  values  of^U)at  ya  and  v\*v  . 

Equations  (9)  and  (10)  are  recursion  relationships  such  that  given 
values  for  viy  and  N4»  which  are  the  values  of  u and  vV  just  above  the 
bottom  we  can  calculate  '-a  and  vA'  at  all  levels  in  the  water. 

Since  we  are  interested  in  normalized  values  of  va  over  the  water 
column  we  can  select  \Ay  to  be  any  arbitrary  value  (\A,M  is  convenient.). 

For  an  arbitrary  value  of  the  horizontal  wave  number  Kr  which  is  restricted 
by  equation  (3c)  we  can  evaluate  by  equation  (3a).  Then  we  can 
determine  u for  all  levels  by  using  equations  (9)  and  (10)  repeatedly, 
if  the  value  of  V.,.  corresponds  to  a mode,  equation  (3b)  will  be 
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satisfied  at  the  surface  of  the  water.  Each  mode  has  at  most  one  such 
solution  for  a given  frequency.  There  is  a low  cut-off  frequency  for  each 
mode  so  that  at  frequencies  below  the  cut-off  frequency  equation  (8b) 
cannot  be  satisfied. 

After  finding  the  amplitude  distribution  of  a mode  it  is  necessary 
to  define  the  mode  number.  For  a finite  frequency  the  mode  number  is 
equal  to  the  number  of  nodes  in  the  amplitude  distribution.  Thus,  the 
first  mode  has  a node  only  at  the  surface.  A representation  of  the 
amplitude  distribution  of  the  first  mode  is  shown  in  Figure  1. 


B.  Ray  Equivalent 

Corresponding  to  the  definition  of  a mode  discussed  above  is  a 
more  physical  approach  in  which  the  ray  equivalent  of  the  solution  to  the 
wave  equation  is  considered.  This  approach  has  been  discussed  previously, 
References  5 and  6,  and  in  Reference  7. 


For  simplicity  let  us  consider  a two  layer  medium  of  constant  water 
depth  H,  density^p,  , and  sound  velocity  cv  , lying  over  a infinite 
bottom  of  density  and  sound  velocity  c*.  as  shown  in  Figure  3.  At 
large  ranges  from  a point  source  we  may  consider  sound  to  be  propagated 
by  plane  waves.  It  is  clear  from  Figure  3 that  for  certain  waves  whose 
direction  is  defined  by  an  angle  Q>  there  will  be  constructive  interfer- 
ence between  it  and  a plane  wave  which  undergoes  one  more  bottom  and 
surface  reflections.  In  order  for  constructive  interference  to  exist  the 
phase  difference  between  points  A and  B,  Figure  3 must  be  10r\-\)TT  degrees, 
or  it  must  satisfy  the  equation: 


an  [ _ 


a. 


Cut© 


*X  G 


) _ 


tt  3 a.C'A-v'jn 


or 

— [ a>rtCou©"\  - t - T»  ’ iU-»)  TT  , . 

X.  \_  -1  (11) 

where  X*  is  the  wavelength  of  the  preferred  mode,  fc  is  the  phase 
change  undergone  by  a plane  weve  upon  bottom  reflection,  r\  is  the  mode 
number,  and  the  phase  change  upon  reflection  from  the  water  surface  is 
assumed  to  be  — v:  . If  the  sound  velocity  in  the  water  layer  varies 
with  depth  the  first  term  of  equation  (11)  would  be  different  from  that 
given  above.  However,  the  discussion  below  applies  to  either  case. 

If  for  a given  wavelength  and  angle  0 there  is  constructive 
interference  between  plane  waves  suffering  different  numbers  of  bottom 
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and  surface  reflections,  then  propagation  consists  of  a series  of  upgoing 
and  downgoing  waves  as  shown  in  Figure  4.  The  left  hand  term  of 
equation  (ll)  is  the  phase  change,  , undergone  in  the  2 direction 

when  a ray  makes  a surface-bottom-surface  cycle.  For  finite  frequencies 

TT  > fc.  * O 

as 

Therefore,  the  phase  change,  & , undergone  in  the  i.  direction  over  the 
water  depth  is  limited  by  the  following 

£».  (12) 

The  pressure  is  zero  at  the  surface,  the  phase  change  upon  reflection  is 
— Tt  , and  the  direction  of  propagation  is  reversed  upon  reflection  from 
the  surface.  Therefore,  the  sound  field  in  the  vertical  direction  is  the 
sum  of  two  sine  waves,  representing  the  upgoing  and  downgoing  waves  in 
the  t direction.  These  waves  are  shown  for  the  first  two  modes  in 
Figure  5.  Because  of  equation  (12)  the  number  of  nodes  in  this  amplitude 
distribution  is  equal  to  the  mode  number.  Thus  there  is  a correspondence 
between  the  definitions  of  mode  number  in  the  solution  of  the  wave 
equation  and  the  ray  equivalent  solution. 

When  the  wave  equation  is  solved  numerically  values  of  are 

obtained.  4U)  is  given  in  equation  (7)  by 

^ 

where  V*-  is  given  in  equation  (5)  by 

Vr  1 ^ »»*vn0 

Therefore  given  positive  j/M  one  can  determine  from  equations  (13) 
and  (14)  the  cosine  of  the  angle  of  inclination  of  the  equivalent  ray 

Cos.  © * ^ C$CV>  (15) 

as  a function  of  1 . 

If  the  ray  between  two  points  and  "2.v  is  continuous  then 
• 1.  - 1 1 may  be  given  in  terms  of  the  horizontal  distance 

and  one  particular  value  of  the  tanQ,  over  the  path. 
This  relationship  is 


(13) 

(14) 


Vc\V\  Q 
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^ ^>v\€)  » c. 

^ i . C- 

c_ 

c - 

C-  -» C-cv>  0 


where  c-*  the  vertex  velocity  = C- 

V w\© 


If  the  value  of  & does  not  vary  appreciably  between  1,  and  '4. 
we  can  approximate  tanO  by 

0 v V b A. 


A 

where  and  are  taken  at  a, 
(Reference  8) 


and  -a*,  respectively.  Thus 


-i_  C-\  V C\ 

^ C-CA^  r 


(16) 


then 


where  C-^  - vertex  velocity 

c.,,Ci=  are  the  sound  velocity  at  ^ and 
^ >U'A=  angle  relative  to  normal  of  ray  at  T,  and  lx. 


Thus  given  values  of  one  can  construct  a ray  equivalent.  If 

vertexing  takes  place  the  depth  1*  at  which  this  occurs  is  the  level 
whose  value  of  sound  velocity  equals  . 


C.  Phase  Velocity  and  Group  Velocity 

The  phase  velocity  V|p  is  given  by  the  relationship 

M t _ 

* (17) 

where  0 is  the  direction  of  propagation  of  a plane  wave  where 
the  sound  velocity  has  a value  C.  Equation  (17)  can  also  be  written 

(i8) 

where  is  the  vertex  velocity 

The  group  velocity  can  be  considered  from  two  points  of  view, 
lrst  group  velocity  >ia  may  be  considered  as  a measure  of  the  speed 
of  propagation  in  the  horizontal  direction  of  a number  of  frequencies 
in  a band  centered  about  ^ may  be  given  by  (Reference  9) 


\J 


S - 


(19) 
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or  by  (Reference  10) 

v 

<—  c_ 

It  can  be  seen  that  this  approach  to  the  calculation  of  group 
velocities  involves  the  calculation  of  derivatives.  This  is  not 
desirable  in  a computer  program  since  this  produces  inaccuracies  and 
makes  it  necessary  to  obtain  an  unnecessarily  large  number  of  values  of 
group  velocity  as  a function  of  frequency. 


(20) 


Tolstoy  (Reference  11 ) has  used  a general  theorem  by  Biot 
(Reference  12)  to  show  the  equivalence  of  ^ in  equations  19  and  20  to 
the  rate  of  energy  transport  in  the  horizontal  direction.  The  group 
velocity  is  given  by 


where 


\ \)_ 

CT 

(21) 

/ I**  . » 

\i  * 

( 

(22) 

cr  - 

(23) 

where  and  c.  are  respectively  the 
at  , and  <$  is  given  in  the  equation 

X.  JL  , - v ><Ct  | 

$ * <=. 


density  and  sound  velocity 


(24) 


where  is  the  displacement  potential 


in  equation  (l). 


Since  by  equation  (2) 

a-t1- 


(2) 


and  u is  the  value  of  pressure,  normalized  to  maximum  amplitude,  as 
a function  of  then  by  equations  (2)  and  (24) 


V-X 


* 


where  is  the  normalized  value  of  $ . 


(25) 


Thus  by  (25)  we  can  obtain  once  >j^  is  known  since  we  are  only 
interested  in  the  normalized  values  of  va  and  for  given  . 


Given  J> . in  the  water  and in  the  bottom  and  u normalized  to 
maximum  pressure  amplitude,  in  order  to  obtain  , the  normalization  of 


9 


where  \ and  v>  signify  water  and  bottom  and  signifies 

unnormalized.  Since  the  maximum  value  of  vA  lies  in  the  water,  since 
i*  greater  than  , and  because  vAt-*)  is  normalized  with  respect 
to  maximum  amplitude  we  multiply  both  expressions  in  (25a)  by  J3,  to 
obtain  normalized  with  respect  to  maximum  amplitude  so  that 

o , (25b) 

JV  P* 

wf  is  the  normalized  value  for  d>  and  may  be  used  in  place  of  <$>  in 
equations  (21-23). 


D.  ExcitationPressure  and  Propagation  Loss 

The  sound  field  produced  by  a simple  harmonic  source  in  a two-layered 
half-space  (shown  in  Figure  6)  with  a free  surface  at  o and  the 
boundary  between  two  fluids  at  l-V  is  given  by  the  solution  of  equation 
(1) . This  is  given  in  Reference  (9)  by 

$ ■ (26) 

where  v~  = horizontal  range 
= radial  frequency 
w\  = mode  number 

= horizontal  wave  numbery^  for  mode  m 


where  Jl  is  the  water  density  at  the  source 

v^«^x)is  the  normalized  displacement  potential,  a function  of 
depth 

is  the  source  depth 
t is  the  receiver  depth 
wv  is  the  mode  number 

\p^  is  the  excitation  function  given  by 

^ j><  c.c  S V V (2S) 

^ >*v\  'Twvvx 


where 


and 
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_A  is  the  water  density 

is  sound  velocity  at  the  source 
5>  is  the  power  output  ofsn  omnidirectional  source1 
is  the  horizontal  wave  number 


^ 1L 


where  ^ is  the  normalized  displacement  potential. 


(29) 


If  the  source  produces  a unit,  sound  pressure  level  then  the  following 
relationship  (Reference  9)  must  be  satisfied 


Ha  ( S<-~ 


K XT?  J 


\ 


(30) 


Substituting  (30)  into  equation  (28)  we  obtain  the  exitation 
pressure  ^ such  that 

p„  ■ (31) 

WN  ( VIva 

The  exitation  pressure  is  the  sound  pressure  amplitude  produced  by 
a source  which  generates  a unit  source  pressure  level  at  unit  range 
when  both  source  and  receiver  are  at  antinodes.  It  is  essentially  a 
measure  of  the  source  level  of  mode  m for  a unit  source. 

From  equation  (26)  we  determine  the  pressure  amplitude  character- 
istics of  the  sound  field  and  this  amplitude,  ifv  , is  given  by 


P--  «*> 

^ it"  J yi* 

Since  ft*  is  the  sound  pressure'amplitude  at  range  r from  a generator 
with  unit  source  level,  the  value  of  propagation  loss  V,r  at  range  V-  for 
a given  depth  of  source  and  receiver  is  given  by 


\-r  ' - lO  \= 


(33) 


It  can  be  seen  from  equations  (27)  and  (32)  that  once  ^is  known 
one  can  easily  determine  the  effect  of  the  source  and  receiver  depth 
on  the  sound  field  at  a given  range.  If  the  source  depth  is  such  that 
» where  1*  is  the  source  depth,  is  a node  the  mode  v*\  will 


LThe  quantity  is  represented  by  the  symbol  of  capital  TV  in  Reference  9. 
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be  suppressed  in  the  sound  field.  Likewise,  if  is  such  that 
is  an  antinode  the  sound  field  of  mode  m will  be  greater  than  at  any  depth 
for  which  is  less  than  ^vi  . The  preceding  also  applies  to  a 

discussion  of  the  effect  of  the  receiver  depth  upon  the  sound  field. 

It  can  also  be  seen  from  equations  (27),  (28),  and  (29)  that  the 
pressure  amplitude  determined  does  not  depend  upon  the  normalization 
of  <$>CW  . 

For  small  attenuation  of  individual  modes  equation  (32)  may  be 
rewritten  to  include  losses  at  the  boundaries  as  a function  of  range 
for  mode  v»\  so  that 

* VL  Vo  V_'*_VVN'C-  - ^ 


>1 


(34) 


where  is  given  in  db  of  loss  per  unit  increment  in  range. 


DESCRIPTION  OF  PROGRAM  SI 441 

Program  SI 441  uses  normal  mode  theory  to  predict  acoustic  propagation 
over  a flat  homogeneous  ocean  bottom.  The  program  may  be  used  with  any 
of  three  options,  each  of  which  provides  different  information  about  the 
sound  field  in  a medium  for  a given  frequency,  velocity  profile,  and  mode 
number. 


The  three  options  provide  ns  output  the  following  calcomp  plots: 

A.  The  first  option  produces  two  plots  for  each  mode  analyzed. 

One  plot  gives  pressure  normalized  to  the  maximum  amplitude  as  a function 
of  depth.  Another  plot  gives  the  ray  equivalent  of  the  mode. 

P.  The  second  option  produces  two  plots  for  each  mode.  One 
plot  gives  three  quantities:  phase  velocity,  group  velocity,  and 
excitation  pressure.  These  quantities  are  plotted  as  a function  of 
frequency.  Another  plot  gives  the  angle  of  incidence  of  energy  at  the 
two  boundaries  for  the  given  modes.  These  angles  are  plotted  as  a function 
of  frequency. 

C.  The  third  option  produces  a plot  of  propagation  loss  versus 
range  for  any  combination  of  modes.  These  plots  can  be  produced  for  any 
source  or  receiver  depth  or  frequency. 
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For  any  of  the  three  options  used,  two  plots  of  the  velocity 
profile,  shown  in  Figures  7 and  8,  are  produced.  The  first  plot 
(Figure  7)  shows  the  sound  velocity  in  both  the  water  and  the  bottom.  The 
second  plot  (Figure  8)  shows  in  more  detail  via  an  expanded  velocity  scale 
the  sound  velocity  in  the  water. 

The  format  of  the  input  data  is  shown  in  Tables  I to  IX.  In  Table  II 
it  is  shown  how  to  obtain  either  of  options  A,  B,  or  C.  The  mechanics  of 
the  individual  options  are  described  below. 


Option  A 

In  option  A the  numerical  solution  to  equation  (4)  or  (6)  is  found 
subject  to  the  boundary  conditions  given  by  equations  (8a),  (8b),  and 
(8c).  Then  the  ray  equivalent  of  this  solution  is  obtained. 

For  a given  velocity  profile  inputed  (Table  IV)  tne  sound  velocity 
is  calculated  at  N (Table  V)  levels  equispacea  between  the  surface  and 
bottom  by  interpolation  between  the  values  given. 

The  value  of  the  horizontal  wave  number  W,.  is  varied  subject  to 
the  restrictions  given  in  equation  (8c)  and  for  each  value  of  Vv-^vj&iis 
calculated  over  the  water  column  of  N equispaced  levels  by  means  of 
equations  (9)  and  (10).  The  values  of is  restricted  such  that 
never  exceeds  (Table  V).  Mode  Vv  must  have  VA  zero  crossings. 

After  is  calculated  for  a given  V*.  > is  incremented  by  OVtA  so  as 
to  obtain  the  smallest  possible  value  of  at  the  surface* 

If  the  conditions  necessary  for  the  existence  of  a given  mode 
cannot  be  met,  the  statement  "No  Mode  Found"  is  printed  out. 

Once  V<D  has  been  found  for  a given  mode  the  ray  equivalent  can  be 
found  by  equation  (16).  Sample  calcomp  plots  of  the  amplitude  distri- 
bution normalized  to  maximum  amplitude  and  ray  equivalent  are  found  in 
Figures  1 and  2. 

For  a large  negative  velocity  gradient  high  frequency  sound  is 
trapped  near  the  oceans  bottom.  Under  these  circumstances,  it  is 
difficult  to  obtain  a good  approximation  to  the  boundary  condition  at 
the  surface.  Then  it  is  necessary  to  decrese  the  value  by  which  V,-  is 
incremented.  This  is ■ done  by  using  a large  value  (up  to  about  10)  of 
IEX  (Table  II).  A large  value  of  IEX  will  increase  the  program  time 
since  it  decreases  the  increment  of  k*.  by  a factor  of  ^q-IEX.  Under 
most  circumstances,  a value  of  IEX  of  zero  is  adequate. 


Option  B 
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In  option  B group  velocity,  phase  velocity,  excitation  pressure, 
and  angle  of  incidence  at  the  boundaries  relative  to  the  normal  to  the 
boundary  are  found  over  any  frequency  range  (Table  VI).  For  a given 
frequency,  group  velocity  is  determined  by  equations  (21),  (22),  and 
(23);  phase  velocity  by  equation  (IS),  and  excitation  pressure  by 
equation  (31).  The  angle  of  incidence  at  the  boundaries  is  given  by 
equation  (15).  If  the  ray  vertexes  before  striking  a boundary,  the 
angle  of  incidence  is  given  as  90°.  Sample  plots  obtained  from  Option  B 
are  shown  in  Figures  9 and  10. 


Option  C 

In  Option  C propagation  loss  for  a source  level  at  a one  yard 
reference  is  obtained  as  a function  of  range  for  any  given  frequency  and 
combination  of  modes.  Propagation  loss  may  be  calculated  by  equations 
(32)  and  (33).  The  ranges  over  which  loss  is  plotted  and  source  and 
receiver  depth  are  inputed  by  a card  described  in  Table  VII.  Values  of 
Ow  are  inputed  as  shown  in  Table  VIII  and  the  modes  which  make  up  the 
sound  field  are  inputted  as  shown  in  Table  IX.  A sample  calcomp  plot 
is  shown  in  Figure  11. 


CONCLUSIONS 

Program  S1441  Is  designed  to  calculate  and  plot  many  quantities  of 
interest  in  the  study  of  a sound  field  in  an  ocean  bounded  by  flat 
parallel  boundaries.  The  problem  of  acoustic  propagation  is  approached 
from  the  standpoint  both  of  physical  and  ray  acoustics. 

It  is  possible  at  high  frequencies  when  the  velocity  profile  has 
more  than  one  "channel"  to  obtain  solutions  of  the  differential  equation 
(6)  which  are  physically  unrealizable.  This  can  only  be  detected  by 
observing  the  values  of  y*'  which  are  3hown  in  the  printout  from  this 
program.  If  two  depths  which  have  positive  values  of  which  indicate 
propagation  of  rays  are  separated  by  negative  values  of  which 

indicate  a "shadow  zone"  one  has  a physically  unrealizable  situation. 
Since  the  velocity  profiles  taken  over  the  BIFI  range  show  in  general 
a monotonically  decreasing  or  increasing  sound  velocity  with  depth  this 
limi table  does  not  severely  handicap  analysis.  Howevpr,  caution  should 
be  taken  in  analysis  when  complex  velocity  gradients  are  used. 


WILLIAM  0.  KANABIS 
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TABLE  I 
Heading  Card 
CARD  1 


INPUT 

FORMAT 

(110,  F10.0,  4110) 

PARAMETER 

COLUMNS 

NUMV 

1-10 

Number  of  velocities  in  profile 

VEL  1 

11-20 

Velocity  at  origin  of  calcomp 
plot  of  velocity  profile 

IV0P 

21-30 

If  IRP  = 0 

When  IVOP  = 0 Option  A 
When  IVOP  = 1 Option  B 


IRP 

31-40 

When  IVOP  = 1 

If  IRP  = 1 Option  C 

IVRP 

41-50 

IVRP  = 1 Velocity  profile  not 

plotted 

IVPP  = 0 Profile  plotted 

IEX 

51-60 

IEX  changes  increment  of  Kr  by  a 
factor  of  ^p-IEX.  Values  from  0-! 

TABLE  II 


CARD  II 
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FORMAT  (5F10.3) 


INPUT 

PARAMETER 

COLUMNS 

ZM 

1-10 

Water  depth  (ft) 

CB 

11-20 

Velocity  of  sound  in  the  bottom 
(ft/sec) 

RO 

21-30 

Density  of  water  (grains  1 cm^) 

RB 

31-40 

Density  of  Bottom 
(grams  1 cm^) 

FSC 

41-50 

Maximum  Depth  Plotted  in  velocity 
profiles  and  plots  in  option  A is 
200* FSC 

TABLE  III 


FORMAT 

(2F10.3) 

INPUT 

PARAMETER 

COLUMNS 

Z(I) 

1-10 

Height  above  bottom  at  which 
sound  velocity  is  C(I) 

C(I) 

11-20 

Velocity  of  sound  in  (ft/sec) 
at  Z(I) 

1=1 NUMV 

in  order  of  increasing  Z 


TABLE  IV 

GROUP  4 
NUMV  Cards 


Li. 


u 
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FORMAT  (110,  2F10.3) 

INPUT 


PARAMETER 

COLUMNS 

N 

1-10 

Number  of  intervals  into  which 
depth  is  to  be  sub-divided  in 
integration  of  differential 
equations 

■ 

UM 

11-20 

Maximum  value  of  vaC-O 

1 

FQ 

21-30 

Frequency  (Hz) 

TABLE  V 

CARD  5 

FORMAT  (4110) 

j 

] 

INPUT 

PARAMETER 

COLUMNS 

1 

ISFQ 

1-10 

Highest  frequency  to  be  analyzed 

IEFQ 

11-20 

Lowest  frequency  to  be  analyzed 

NMOD 

21-30 

Number  of  modes  analyzed 

Mode  numbers  analyzed  =1,  ...  NMOD 

INCF 

31-40 

Decrement  in  frequency  from 

ISFQ  to  IEFO 

Options  A and  C ISFQ  = IEFQ  = FQ 

Option  B Range  of  frequencies  is  selected  by 

ISFQ,  IEFQ,  INCF 


TABLE  VI 
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T"  | - — - I ■ ■ 


FORMAT  (4110,  3F10.3) 


INPUT 

PARAMETER 

COLUMNS 

IRST 

1-10 

Lowest  range  in  feet  at  which 
propagation  loss  versus  range 
will  be  plotted 

IREN 

11-20 

Largest  range  in  feet  at  which 
propagation  loss  will  be  plotted 

IRIC 

21-30 

Increment  in  range  to  be  plotted 
in  feet 

NFS 

31-40 

Number  of  propagation  loss  versus 
range  plots 

ZS 

41-50 

Source  depth  in  feet 

ZRC 

51-60 

Receiver  in  depth  in  feet 

FMI 

61-70 

Increment  of  range  on  the  calcomp 
plot  per  division  in  nautical  miles 

TABLE  VII 


CARD  7 

For  Option  C 


FORMAT  (4F10.5) 


INPUT 

PARAMETER 

COLUMNS 

DD(J) 

1-10 

11-20 

21-30 

DD(J)  is  the  loss  at  boundaries  as 
a function  of  range  in  db  per  foot 
for  the  mode  J. 


J - 1,  NMOD  where  NMOD  is  the 
number  of  modes  analyzed 


NMS  1-10  Number  of  modes  to  be  summed  in 

the  plot  of  propagation  less 
versus  range (NM4NM0D) 


Card(s)  b. 


FORMAT  6110 


INPUT 

PARAMETER  COLUMNS 


MDS(J) 


11-20 


21-30 

31-40 


J = 1,  NMS 

The  value  of  MDS(J)  are  the 
mode  numbers  of  the  modes  to  be 
summed  in  the  propagation  loss 
versus  range  plots. 

MDS(j)^  NM0D 


41-50 


51-60 


Number  of  card(s)  b.  is  the  integer  greater  than  or  equal  to  NMS 

6 

Cards  a and  b form  a set.  There  is  a set  of  cards  for  each  plot  in 
Option  C.  (NPS  sets) 


TABLE  IX 

GROUP  9 
For  Option  C 


Card  a. 


The  number  of  cards  in  this  group  is  the  integer  greater  than  or  equal 
to  the  quantity  NMOD 
4 


INPUT 

PARAMETER 


FORMAT  110 
COLUMNS 
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TABLE  VIII 


GROUP  8 
For  Option  C 
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L A. 

S \ H ^ V 

DOUBLE  PRECISION  F * U . UP 

OOUbLE  PRtClSION  K2,  N2.P2. h»2P»  CZ.G.HZ.DEN 

DOUBLE  PRECISION  S.K20.UK2.6AM.BA 
DOUBLE  PRtClSION  E,V 
oouble  precision  UMAX 
UOuHLt  PREulSION  COS2  • ccc  » cv 

DIMENSION  hED<20) »Z(lO0) .C( l00) »F(50o> » jP(bOo>  .0(500)  . 

1  ZMPUOOO).  XMP(loOO).  C 062(500) . CCC(bOO) 

1 . ZZP<50u>  » 2PP(500)>  VV(500) . DATA(1024)  » UU<500) 

1 »PM(6000)  .DOT(IOO) 

1 » ZR(50U)  * CR(boO)  tKPUOO)  *ZBT(5)  . ZBF(5) 

lrUb(lOO) .PhT( 1000) .PCT(iOOO) .PST < iOQu > . MOS ( 100 ) ,USR| 100). URR(IOO), 
1 PMM(IOO) .FMR(IOO)  .FDS(IOO) 

1 »uVEL(40uo) » PVEL(40u0)»  FQP(4000)»  T 1 (4000 ) . T2(4000) 

COMMON  MS. jM.U 
COMMON  jms 

call  plots  ( data<u,  1024,  o ) 

4 JPM  = 25 

«EAD(3,1)HeD 

1 FORMAT  (20a4) 

READ  (3.112)  NUMV.  VEL1  » IVOP  . IKP  * IVPl  » IEX 

112  FORMAT  ( H0*F10. 0,4110) 

NUM  = 0 

READ (5. 2)  ZM#CB,RO,RB  »FSC 

2 FORMAT  (5F10.3) 

ZMo  = (200.0*FSc  -ZM)/(20.0*FSC) 

ZbT(l)  = ZM 
ZBT (2)  = ZM 
ZbT (3)  = 200.0  *FSC 
ZBT ( 4 ) = -20.0*FSC 
ZbMl)  = o.O 
ZBF (2)  = lu.O 
ZdF(31  s o.O 
ZbM4)  = i.O 
00  3 1=1.100 
Kl)  = 0. 

3 C(i)  = 0. 

1=1 

5 READ  (3.30)  Z(l).ClI) 

30  FORMAT (2Flo.3) 

IF(DABS(Z(l)-ZM)-.0l)o.b,7 
7 I = !♦! 

60  TO  5 

6 1 = 1 

IF  (IVPL.EQ.l)  60  TO  6 
DO  130  K = 1.  NuMV 

CK(NUMV  -K  *1)  = C (K ) 

130  ZR (NUMV  "K  +1 ) = ZOO 

ZR  (NUMV  *1)  = ZM  -0.01 
ZH (NUMV  +£)  = ZM  ♦sO.U 
CR (NUMV  +i)  = C6 

CH  (NUMV  ♦^i)  = Cb 

DO  1B5  NlNY  = 1.  NUMV 
165  ZR (NlNY ) = ZM  - 2R(NINY) 

ZR (NUMV  +3)  = 200.0*FSC 
ZR (NUMV  *4)  = -20,0*FSC 
CR (NUMV  ♦J)  = VtLl 

CH (NUMV  ♦*)  = 50.0 


A-l 
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I 


I 


1 

I 

it 

« 


call  plot  (0.o*o.o,-3) 

call  LINE.  (CR»ZR»NuMV  *2,  1,0,0  ) 

CAcL  LINE  (ZOF»ZBT,2»1,0,0) 

CALL  SYMBOL  1 1 0 , 2b , ZMo » 0 . 1 4 , 6hboT T OM , 0 , 0 , 6 ) 

CALL  Axis  (0.0*0.0,14HSOUNU  VELOCITY, 14, 10. »0.U,CR(NUMV  *3)  » 
1CK (NUMv+4 ) , 10 , 0 ) 

CALL  AXIS  (U, 0,0,0, llHWATtK  DEPTH, 11 , 10. 0,90. 0*ZR(NuMV  a3>  , 

1ZK(NUMV>4) , 10,0  ) 

CALL  PLOT  (lb. 0,0. 0,-3) 

ZH(NUMy  +1)  = iOO,0*PSC 

ZR(NuMy  ♦<>>  s -20.o*PSC 

CRtNUMV  *1>  =VtLl 

CR(NUHv  *2)  510,0 

CALL  LINE  (CH,ZR,NuMv, 1,0,0) 

CALL  LINE  (Z»F,ZBT,2,1,0,0) 

Call  symbol  (10.25,ZMB»0.14.6HBOTTOM,0.0»6  ) 

CAiL  AXIS  (0«0,0.0,14HSOUNO  VELOCITY , 14, 10. »0.0,CR (NUMV  tl)  , 
1CR(NUMV*2) ,10,0  ) 

CALL  AXIS  (0.0,0,0,11H4ATEK  DEPTH, 11, 10.0, 90. 0»ZR(NuMV  +1)  , 

1ZK(NUMV42) ,10,0  ) 

CALL  PLOT  ( lb. 0,0.0,  -3) 

B JP  5 3 

lt)B  *K1TE<4,9)hCD 

9  FOhMAT  (1H1,9X,20A4) 

*RItE<4,10) 

10  FORMAT UHQ,1SX»4HZ, ft, 5X, BMC, FT/SEC) 

11  WRITE(4,12)ZRU)»CR(I) 

12  FOrMAT (1H0,9X,F10»1,Fi2»1) 

IF (DABS (Z ( I )-ZM)-t 01) 14,14,13 

13  I = in 
JP  = JP41 
IF(JP-JPM) il,8,8 

14  RLAD(3,lb)N,  OM  ,FQ 
lb  FOkVAT(I1o,2F1C,3) 

NNN  5 N 

lbl  READ  (3,lbu)  XSFO,  IEFO,  NMoO,  1NCF 
ISO  FORMAT  (4U0) 

IF  (iRP.Nt.l)  GO  TO  120 

KlaD  (3,121)  IRST,  IREN,  IRlC,  t,PS»  ZS,  ZRC  ,FMI 

121  FOkMAT  (4U0*  3F10.3) 

ZS  = ZM  - 2S 

ZHC  5 ZM  - ZRC 

READ  (3,12c)  (OD( j) , Js  1,nMOo) 

122  format  (4kio.S) 

120  DO  1S4  Ms  l.NMOU 

00  lbb  Ls  ISFU,  IeFQ,-InCF 
IF  (IVOP.Ut.D  Go  TO  162 
FQ  s L 

loc  P2  = 6,283id5306 

42  s (P2*Fa)*<P2*Fo) 

N S NNN 
FMoD  5 V 
IRlL  5 o 
IF (N)29.4,c9 

29  K20  = ((42*(CB-C(l))*tCB4C(i)))/(L(l)*C(l)*Co*LlO) 

0K2=  ,U99*n20 
K2  = “N20 
31  X 2 = *2  ♦ jK2 

A-.' 


V 


I 

I 
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It*  (K2~K20) 16*32,32 

32  AhlTEt4,33) 

33  FOhMAT ( 1H0 # liHuO  Hoot  FOUND) 

ItFQ  = u ♦ INCH 

IF  (lVOP.fc.u(.l)  00  TO  lbO 

00  TO  Hi 

lo  FLu  = FtoAT(N) 

H = Z*/  ObLEtPLB) 

ZP  = 0, 

0=1 

1 = 1 

17  IF  (ZP-  Hi+D)  18,13,19 

19  I = !♦! 

Oo  TO  17 

13  CZ  = (C(I)*iZ(UD-ZP)*C(I+l)*(ZP-Z(I)))/(ZtUl)-Z(D) 

CCC<0)  = c* 

0 =-(w2*(c«l)-CZ)*(C(l)*CZ) )/(CZ*LZ*L<l)*C(l> ) 

F ( j)  = K2  - 0 
IMZP-ZM)  <>0*21,21 

20  J=o*l 

FLA  = Fl0AT(0  -1) 

ZP  = ZM*jBLE(FLA)/DBLE(FLB> 

00  TO  17 

21  UID  = 1,0 

UP(1)  =(Ko/Rl3)*0S0RT  ( <W2*(CB-C(1>  )*(Cb+Ctl) ) > / (C  ( 1 ) • <C  < 1 UCB*CB> ) 
1-K2) 

J = 1 

22  OP  = 3 

2b  fla  = float (o  -u 

FLo  = FlOAT(N) 

ZP  = ZM*oBLE(FLA)/DBLE(FLB> 

IF{0-N-1)27,34,27 

27  IF(DA8S(U(J> >-UM)2o*34,34 

28  0 = 0*1 
0P=  0P-*1 
H2=  h*H 

OEN  = 1.0  ♦ (H2*F ( j) )/b,0 

U(o)=  ( (l.u-tH2/3.0)*F(O-l)J*U(O-l)^H*UptO-l) )/DEN 
0P(0)  = Ul.U-'H2*Flu)  ) /J,0 ) *UP (0-1  >-VFtO-mFlO)-lH2*FlO)*F  CO-1)  >/ 
16,o)*H*U< j-l>*  0.b)/0fN 
IF  <0P-0PM)2b,22,2b 
1H  OM  : 0 

call  count 

19b  IF  (M5-P)3l,35,3b 

35  K2  s K2  - oK2 
DK2=  0K2  / 10,0 

IF  (Dk2-  ,Qo0001*K20  *10** (- [EX ) ) 3b, 3b* 31 

36  J = 2 

S = H*0(l)*U(l)/2,j 

52=  H*U(l)*U(l)/(2.0*CCC(l)*CCt(l)  ) 

37  IF (O-OMS) jo,39,38 

38  S = H*U(0)»U(0US 

S2=  H*UlO)*U(U)/(CtC(o)*CCC(0) ) ♦bi 
0 = 0*1 
GO  TO  37 

39  S =(h*U(0)*U(0)/2,0)  ♦ S 

52=  H*U(0)*U(0)/(2,0*lCC(0)*CCC(0)  ) ♦3«! 

IF  (IVOF'.Eo.l)  00  TO  104 

A-» 
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J = 1 
*4  UP  = J 

WHlTE(4,2b)HEU.N.K2»UN.F<J»M 

2b  FORMAT  ( 1H1 , bX , 20A4 , / , ''A . 3HN  =,Ib.2X.4HK2  = .F  lb.  6 . 2X . 4HUM  =»F6.1# 
12a.6HFrEQ  = .Fo.l.^X.bHM  = .l5) 

WHlTE(4,4  f)  ZM.cB.mO.HB 

47  F0kMAT(3X,i,mZMmX  = .F10.3.4HCB  = »F10.b,4HHO  = .F15.6.4HKB  = . 
lFlb.tJ) 
wKlTE(4,24) 

24  FO*MAT(lMo,i7X.4HZ,FT,10X»lHL',loX,5HDU/L>Z.liA.4HF(Z)  ) 

4b  ZP  = ZM4(  FLOAT (J-!)/  FLOAT(N) ) 

Zti  = ZM  - IP 

MK|TE(4»2t>)Z0»0(0)  ,UP<J> .F<j) 

1FIJ-N-I)4b,44.45 

4b  IF (DAUb(Uio) )“04)4b*44»44 

40  J S J41 
JP  = JP+1 

IF  ( JP”JPM)t*b»42»43 
44  WHITE (4*40 ) S 

40  FOkMAT  (Ihu.lOX.  3nS  = .Flb.S) 
lo4  UMAX  =0,0 

NM  S N +1 
DO  104  K=  1*  NM 

IF  ( DABS  (U {K ) ) . UMAX)  104. 104. 103 
10b  UMAX  = DaBS(U(K)  ) 

104  CONTINUE 

00  10b  K=  1 * NM 

10b  UU (K ) = U(K)/UMAX 
00  106  K=  1*  NM 

lOo  Zip (X)  = 2M*(FL0aT(K  -1) /FLOAT (N)  ) 

gam  = DSOHT(K20  - k2 ) 

BA  = (R0*C(1) )/ (2,0*RB*CB*GAM*b  ) 

WHITE (4,41)  BA i GAM 

41  FOhMAT(1Hu,10X»5MB/A  =,  »F15,5.bX. 7HGAMMA  = *Flb,5  ) 

Sb=  R0*S/(UMAX*UMaX)4R0*  R0*U ( 1 ) *U ( 1 ) / < UMAX*UMAX*2 . 0*GAM*RB ) 

SSb=  R0*S2/(UMAX*UMAX)  4 R0*R0*U ( 1 ) «U ( 1 ) / ( UMAX*UMAX 

l*2.0*GAM«Rd*CB«CB) 

IF(IVOP.Eu.l)  GO  TO  lb3 
E i DEXP(-„AM*h) 

ZP  = 0,0 
V = 1.0 
NN  = 1 

ZPP(l)  = 0.0 

Well  = UL(1> 

WHITE  (4.2b)HED.N»K2»UM.F<l»M 
WR1TE(4.47)  ZM.CB.RO.RB 
WRITE (4.46) 

46  FOkMAT(1HO,17X.4HZ»FT,10X.1hU  ) 

49  WHlTE(4.2o)ZP*V 

26  FORMAT  (1H0.F22. 2. FU.  3.F13.,x»Flt>.5) 

50  V = V*E 

NN  = NN  *.1 
VV(NN)  s VV (NN  -l)*t 
ZPP(NN)  = ZP  - h 

ZP  = ZP-H 

if  (Zpp(nn)  4 bo.oi  no.  no.  49 

110  call  plot  < 5,0.0. j.  -b> 

DO  164  N1NX  - 1 , NN 
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184  (MlUX » r 2*  - ,.FP,NINX) 

MNV  = ♦ 1 

DO  183  NlN*  z It  N1NV 
183  2Zp(,\llNM)  - /to  - 22P ( NINw ) 

Vv(NN  M)  z -1,0 
VV^NN  +£)  - u,5 

2PP(\M  *1)  = 200,0*FSt 
2PK(NN  *2)  z -*0,0*FbC 
22H(N  *2)  z 200 , 0*FSC 
22h(N  *3)  z -cO.OeFSc 
UU(N*2>  z -1.0 
UL'(N*3)  = 0.5 

CALL  rtXlb  (0.u,0,U»  ^HAMPLl TUDl  » V»4.0»0.0,liO(N42) »UU(NO) *10.0) 

CAUL  A X 1 b (0.0. 0.0, lbheATLR  OEP1H  IF T) , lb, 10.0,90,0.Z2P(N*2> # 

1 2*P(N  *3)  .10.0) 

call  LINE  (UU,  22k.  N+l.  1*0,0) 

CALL  LINE  (VV,  2PP,  NN,  1.0,0  ) 

CAlL  LINE  (ZOC.ZDT.2,1,0,0) 

CAWL  Symbol  ( 10 , 25 , ZMO , 0 , 14 , bHbOTTOK » 0 . 0 , b ) 
call  Symbol (1.3,9, 50, o , 14»9hamplI TuCE, o , o , 9) 

CALL  SYMBOL (1.7,9. 25,0. 10, 8HVEKbUS,0.0,o) 

CALL  SYVBuL(l.«>,9.o0.0.14,bHDEPTH,0.0,5) 

CAlL  SYMPll(0.o,8.75,0.14,19HFkLQUENCY  HZ. 0.0, 19) 

CAuL  NOMUtH (2. 2, 8, 75,0. 14. FQ, 0.0, -1) 

CALL  SYMBOL ( 1. 8. 50,0. 14. 4HMOOE, 0.0,4) 

CAlL  NUMBlb (2,l»ti,j0.0, 14 ,FM0D, 0 ,0,-1) 

NU.m  = N(J,-i  41 
153  XMB(l)  z 0.0 
Z*V(1>  = 0.0 

NM  = N ♦ l 
DO  170  K=  1 , NN. 

IF  'F(l)  .PT.9,0)  (,0  TO  179 

IF(F(K).8t.0.0)  60  TO  17l 

170  1KTB  = K 

171  C0b2(I«TB  +l)  = F(lRT8  ♦l)*cCC(lRTa  4l ) *CCC < IHTB  4l )/*2 
CV  = CCC(I«Tt>  4l)/DbOHT(1.0  - C0S2(lHTtf  4l)) 

2MK(2)  = m*(CV  - CcCUKTB  >1 > > / (CCC ( IKTh) -CCC ( 1RTO  ♦!)) 

XMPU)  = (2MP(2)  - Zn'.P  ( 1 ) ) * (CCC  ( IHTu  ♦!)  4 CV)/(CV«(  DSoH  T ( C0S2 

1 URTB  +1) ) ) ) 

2M|J(1)  = KlOAT(1RTo)4M  - 2K.P  ( 2 ) 

2M)>(2)  = 2mP(  1 ) 4 2MP(2) 

NM  = NM  -1«TB  41 
DO  172  A = 3 , i.M 
TH1  = 90.0 
KO  = IATD  -1  4K 

C0b2(K0-lJ  z F(kO  -l)*CCC(KO  -l)*CCC(KO  -l)/»2 

COb2(KO)=K(KO)*CCC(KO)*CCC(KO)/e2 

TM2  = ( ACO5 ( SOt'T  ( CoS2  1 NM  41RTU  -1)  > ) >*57.2958 

IMCCC(K0).6E.v.V)  60  TO  142 

XMKK)  = Xv.P ( K -1)  4H*  (CCC  (KO  -1>  ♦ CCC(KO)  )/(LY*(OSQKT(COS2(aO-1 
1)  ) ♦ Ub(JWT(C0S2(K0)  ) ) ) 

172  2mh(a)  z (AO  -1 ) *H 
N z N.M  -l 

40  TO  173 

179  COb?(l)  = ABS(Hl)  )*CCC(1)*CCCC1)/a2 
CV  z CCC ( 1 > /CSORT ( 1 . 0 - COSZ(l)  ) 

DO  140  Kz  2,NM 

COa2(K  -1)  = AL'S  ( 4 i K -1))4CvC(a  -1)*CCC(A  — 1 ) / A2 


NUSL  Tech  Mano 
2211-296-69 

COi,2(K)  = ABb  (F  (K ) ) *CCC  (K  ) *CCC  (K ) /W2 
TMl  =(AC0b(bORl  (C0b2(l) ) ) W57.2958 
TH*  =UCOb(bOST  (COi2(WM)  ) > )*57.2958 
IF  ( CCC(K)  .&E..U)  60  TO  14  2 

XMP(K)  = amP<K  -1)  ♦ m*(CCC(K-1)  ♦ CCC(K))/<CY*(DSQKT(C0S2(K-1))+ 
1 USGRT(C0b2<*)  > ) ) 

140  ZMP(K)  = (K-1)*H 

175  00  141  K=  l.N 

XMHJNM  *K)  = 2 • 0*XmP  (NM)  - XMP(NM  - K) 

141  ZMP(NM  +K)  = ZMP'NM  -K) 

60  TO  143 

142  ZMp(K)  = ZmP<K  -1)  ♦ H*  (C  V - CCC (K-l ) ) / (CCC <K)  - CCC(K  -1)) 

Th2  = 90.0 

XMP(K)  = Xmd(K  -1)  ♦ iZMP(K)  - ZMP(K-l) )*(CCC(K-1)  ♦ CV)/(CV*( 

1 OSORT(COb2<K  -1) ) ) ) 

N = K -1 
NM  = K 
DO  144  K=  l.u 

XMP(N+1+K)  = 2 • 0*XhP (N+l ) - XMP(N+1-K) 

144  ZMP(NM  ♦ *)  = ZMP(NM  -K) 

143  XMP(2*NM)  = 0,0 

182  XMp (2*NM  +1)  = 1000. 0 

NIZZ  = 2*Nm  -1 
DO  186  Mu*  = 1,  NlZZ 
lttb  ZMP(NINZ)  = ZM  - ZMP(NINZ) 

ZMp(2*NM>  = 200 • 0*FSC 

ZMP  (2*NV.  +1)  = -20.0*FSC 

WRITE  (4,145)  TH1 , TH2 
IF  (IVOP.tO.l)  60  TO  157 
14b  FORMAT  (2F10.3) 

CALL  PLOT  (l5.O,0,u,-3) 

CALL  LINE  (XMP, ZMP, NM  +N,l,o»0) 

CALL  LINE  (ZBF.ZBT, 2*1,0, 0) 

CAtL  SYMBOL  (10.25,ZMB,0.l4,6HbOTTOM,U.O*6  ) 

CALL  AXIS  (0.0,0.0,10HRAN6E  (FT) , 10 , 12. 0, 0 . 0 , XMP(2*nM)  , 
1XMP(2*NM  ♦i),10.0) 

CALL  AXIS  ( 0 • 0 • 0 , 0 , 16H  OEPTH  (FT) ,1b, 10. 0,90.0. ZMP(2*NM) , 

1 ZMP (2*NM  +1) ,10.0  ) 

CALL  SY»*Bol(5.°»9.sO,U.14,14HRAY  EGO1VAlENT»0.0,14) 

Call  SYMl)oL(5.0,9.25,0.l4»19HFF<tQUENCY  rlZ,0.0,19) 

CALL  NUmBEk (6.4,9.2b,0.14,FQ»0.0,"l) 

CALL  SYMBOl<5.S,9.o0»0.14»4HMODE»0.0»4) 

CAlL  NUMBER (6.2,9,00,0.14, FMOD , 0 . 0 , -1 ) 

CALL  PLOT  (15.0,0.0, -3) 

IF  (IVOP.Nt.i)  Go  TO  154 
157  T1(L)  = TH1 
T2(L)  = TH2 

FMm  = SORT  (P2*F'Q  / (30 ,48*CV ) ) 

PM(L)  = Si)rT(P2  ) * RO/  (FMH*2 . 0*  (SS)  ) 

6VtL(L)  = sS/(LV*SbS) 

FOP (L)  = L 

•RITE  ( 4,166)  PM(L)  , FMh,  SS  , 6VtL(L) 

16o  FOkMAT(4F1o.3) 

IF  (IRP.EO.1)  PM(M)  = PM(L) 

FMR(M)  = P2*F0/CV 
IZi  = ZS*Fl8/ZM  +1.0 

IZK  =ZRC*FlB/ZM  >1.0 

USk(M)  = U(IZS)/UM„X 
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URR(M)  = U(I*R)/UMaX 
PMm(M)  = P,y ( M ) * USR ( M ) *URR ( M ) /RO 
155  PVtL(L)  = O'* 

If-  (IRP.EO.1)  <*0  To  154 

160  IDF  = (ISFQ  “ IEFG)/INCF  +1 
LIF«  = IDF  -1 

00  161  LLL=  1*  IDF 
ILL  = (LLL  -1)*INCF 
FQp(LLL)  = FwP  ( IeFQ+ILL) 

GVeL (LlL)  = 6vEL(IEFQ+ILL> 

Tl(LLL)  = Tl  ( IEFq+IlL) 

T2(LLL)  = T2(IEF*+IlL) 

PM  (LLL)  = pM(  It.FQ  ♦ ILL) 

161  PVEL(LLL)  = PVEHIEFQ+ILL) 

*IKlTE(4rl6j)  (PV El(J)»  GVEL<J),  FQP(j),  PM(J),  JS  i,IDF) 

165  FOkMAT  (4F10.5) 

FQp(LlFQ  +2)  = 0.0 

FG>P(LIF0  +3)  = ISFQ/lO 
PVELCLIFQ  +2)  = 5500.0 

PVEHLIFQ  +3)  = 250.0 

GVtL(LIFQ  +2)  = 5500.0 

GV£L(LIFQ  +3)  = 250.0 

PM(LIFQ  +2)  = 0.0 

PM(LIFQ  +3)  = 0.1 

CALL  PLOT  (15.0,0.0, -3) 

CAlL  LINE  (FQP,PVEl,LIFQ  +1,1,10,14) 

CALL  LINE  (FQP*GVEl»LIFQ  41.1.10*26) 

CALL  LINE  (FQP*PM,lIFQ  +1.1,10.4) 

CALL  AXIS  (0.0.0.0.14HFREGUENCY  (H2) ,14,10.0,0.0»FQp(LIFQ  +2). 

1 FqP(LIFQ  +3) >10.0) 

CAlL  AXIS  (0*0. 0.0, 27hSOUND  VELOCITY  (FT  PER  SEC  1.27.10. 0*90.0. 

1 PVEL  (LIFQ  +2).PvEL(LIFQ  +3),10.0) 

CAlL  AXlS(iO. 0.0.0, 18HEXITATION  PRESSURE. -18, 10. ,90.0»PM(LIFQ  +2). 
1 PM(LIF0  +3) »10,0) 

CAlL  SYMBOi_(2.4,9,50.0.14»38HSOUND  VELOCITY  AND  EXCITATION  PRESSUR 
IE, 0.0. 38) 

CALL  SYMBOL (4.7,9.25,0.10. 6H VERSUS, 0.0,6) 

CAlL  SYMBOl(<*.4»9.O0»O.14,9hFREGUENCY,O.0,9) 

CALL  SYMBOL(‘f.5,8.75,0.14,4HMODE»0.0,4) 

CAlL  NUMBEk(5.2,8.75,0.14»F*Od,0.0,-1) 

CAlL  SYMBOL (4.7,8.50,0.10, 14HPHASE  VELOCITY.0.0,14) 

CAlL  SYMBOL (4»2» 8.50, 0.10, 14 »0,0»-l) 

CAlL  S YMBOL ^ ,8.25,0.10, 14HGR0UP  VtLOClTY»0.0,14) 

CAlL  SYMBOL (A. 2, 8. 25, 0.10, 28,0. 0,-1) 

CALL  SYMBOL (4.7,8.00,0.10, 19HEXC I T AT ION  PRESSURE. 0.0, 19) 

CAlL  SYMBOL (‘♦.2, 8. 00,0.10, 4,0. 0,-1) 

TKLIFQ  +2)  = 9o,0 

TKLIFQ  +3)  = -lo.O 

T2(LIFG  +2)  = 90.0 

T2(LIFQ  +3)  = -lO.o 

CALL  PLOT  (15.0,0. 0,-3) 

CALL  LINE  (FOP,  T1,LIFG  +1,1»10,4) 

CAlL  LINE  (FOP,  T2,LIFG  +1,1*10,5) 

CALL  AXIS  (0.0.0.0.14HFREGUENCY  (HZ) ,-14,10. ,0.0, Fup(LlFQ  +2), 

1 FoP (LIFQ  *3), 1.U0) 

CAlL  AXIS  ( 0 . 0 , 0 « 0 , 28HANGLE  OF  INCIDENCE  (DEGREES) ,26,  9.0,90.0, 

1 TKLIFQ  42),  TKlIFu  +3),l0.0) 

CAeL  SvMool (3.7,9.30,0.14, 18HANGLE  OF  INCIDENCE, 0.0, 16) 
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Call  SYMBOL <4. 7,9. 25,0. 10*6HVERSUS,U.0*o> 

CALL  SY*,Bul<‘*.4*9.00»0.14»9hFREUUENCY*0.0,9) 

CALL  SYF'UOl (4, 5*8.75,0.  14*4HM0UE #0,0*4) 

CALL  NUMBe.w<5,2,8.  75,0. 14, FmOO,0. 0,-1) 

CAlL  SY*'BOu(‘*.b»b,50*0.10,25HBQTTOM  ANGLE  OF  INC1DEnCE*0.0«25) 

CALL  S YMBOu  (4, 0,8, 50*0.10*4, 0,0*"1) 

Call  Symbol (4, 5, 8, 2b, 0.10* 2&HSUHFACL  ANGLE  OF  INCIDENCE, 0.0, 26) 
CALL  SYMBOL {4«0*8«25»0. 10*5,0, 0,-1) 

154  CALL  PLOT  (l5.Q,0.o»-3) 

IF  (IRP.Nt.l)  GO  TO  111 
123  UO  124  IP-  1,NPS 
READ  (3,125)  NMS 

125  FORMAT  (Hu) 

REaO  (3, 12b)  (MOS(J),  J=1.NMS) 

12b  FORMAT  (6110) 

00  131  IR  = IRST , IREN*  I«IC 
IREO  = (IR  -IRST)/lRIc  +1 
PCT(IREO)  = 0.0 

131  PST (IREO)  s 0.0 

DO  127  Ig=  1 , NMS 
Ms  MOS(IQ) 

00  129  IR  = IRST*  IREN*  IRIC 
RR  s IR 

IREO  = (IR  -IRSD/lRIC  +1 
OOT(M)  = lu.0**(-DU(M)*RR/20.0) 

PST (IREO)  = PST ( iREO)+PMM(M)*UDT (M)  *SIN(FMr(M)*RR-P2/B.0> 

PCT(IRE0)=  PCT(1RE0)+PMM(M)*DDT(M)  *COS(FMR(M)*RR-P2/6.0) 

129  CONTINUE 

127  CONTINUE 

DO  132  IR  s IRST*  IREN*  IRIC 
RR  = IR 

IREO  = (IR  -IRSTJ/lRIC  +1 

PPT(IREO)  = SQRT(PCT(IREO)*PCT(IREO)  ♦ PST ( IREO) *PST( IREO) )*RO/ 
1SQRT (RR/3, 0 ) 

PPT(IREO  ) =-20.0*ALOGlO(PPT(lREO>) 

132  RP( IREO  ) = FLOAT (IR)/6000.0 

WRITE  (4*128)  (PPT(J),  J=  1 * IREO  ) 

128  FORMAT  (10F10.5) 

RP(IRE 0 +1)  = 0,0 

RP( IREO  +2)  = FMi 

PPKIREO  +1)  = 3u.O 

PPT  < IREO  +2)  s 10. 0 
2S  s ZM  - IS 

ZRC  = ZM  - ZRC 

CALL  LINE  (RP»PPT*IRE0.1»10.4> 

CALL  AX IS(o, 0*0.0* 13HRANGE  (MILES) ,13*20. 0,0.0, HP( IREO  +1 ) * 

1 RPdREO  +2)  *10,0) 

CALL  AXIS(o.0*0,0*21HPROPAGaTION  LOSS  (Ub) *21*10.0* 90.0, 

1 PPT ( IREO  ♦1>*PPT(IRE0  ♦2)»l0.0) 

CALL  SYMBOL(0.9.9.bO*0,14*16HPROPAGAT10N  LOSS, 0.0, 16) 

CALL  SYMBOL (9.7,9.25*0. 10 *6HVERSUS, 0.0 *6) 

CALL  SYMBOL (9.7,9.00 *0. 14* 5HRANGE* 0.0 *5) 

CALL  SYMBOL(8«9*8.75*0.14*lBHFRt.OUeNCY  HZ. 0.0. 18) 

CALL  NUMBER (10,3?8,75».14*FQ*0,0.-1) 

CAlL  SYMBOL  (8.b,6.b0»U.14*20NSOURCE  UtPTH  FT. 0.0, 20) 

CALL  NUMBER (10. 5*8,50,0. 14, ZS,0.0*-1) 

CAlL  SVMBOL(8*t>*8,25,0.14,22HRECEIVLR  UtPTh 
CAWL  NUMBER ( 10. 6*8, 25, 0.14, ZRC, 0. 0,-1) 
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CAUL  SYMBOL  (8«6*8,0»0.l4* 5Hm00LS 1 0 . 0 » b ) 

00  187  Js  1»NMS 
POi(J)  s MOS(J) 

187  CALL  NOMB&w  (9.1  ♦0,S*J.8.0»0.19#Fl/S(J)»0.0.-l) 

call  plot  us.o»o.o#-j> 

124  CONTINUE 

call  plot  <i5.o,o.o»-3> 

111  CALL  PLOT  cl5.Q»0.0#-3> 

End 
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SUbHOUTINE  COUNT 
DIMENSION  u(b OC) 
DOUBLE  PMtClSION  u 
COmMON  MS»JM*U 
common  JMS 
MBsO 
0=1 
IS=1 
JMS=1 

b IF  (U(J))  1#*#3 

1 ISisiS 
ISsl 

7 IF  (IS-IS1)  *»,2,4 
4 MS=MS*1 
JMS=J 

2 J=u+1 

IF  (J-JM-1)  5#6.5 
1 ISl=lS 
ISsO 
60  TO  7 
b RETURN 
END 
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